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Molecular-Dynamics Simulation of Vulcanian Eruption 

Satoshi Yukawa* and Nobuyasu Ito^ 

Department of Applied Physics, School of Engineering, The University of Tokyo, 
7-3-1, Hongo, Bunkyo-ku, 113-8656, JAPAN 

Vulcanian explosive eruption, which is a nonlinear and nonequilibrium abrupt dynamics of 
magma-gas mixture, is modeled by a two-component Lennard- Jones particle system. Molecular- 
dynamics simulation of a shock-tube experiment gives consistent results with a explosive erup- 
tion picture of volcanology; Shock wave and expansion wave are reproduced. In addition bubble 
nucleation of a gas component in the magma melt and spinodal-like decomposition are observed 
in the simulation. The result is also compared with a continuum hydrodynamic model; Quali- 
tative features of continuum dynamics are reproduced by the present model. We find that the 
particle description of dynamics is an effective method in such kind of abrupt dynamics. 

KEYWORDS: Vulcanian eruption, Lennard- Jones particle, molecular dynamics simulation, shock tube 



Volcanic eruption is complicated physical phenomena 
and the physical understanding has not been well estab- 
lished yet; The problem is to understand nonlinear and 
nonequilibrium dynamics of magma-gas mixture accom- 
panied by phase transitions.^"^ Existence of gas, which is 
mainly H2O, is sometimes forgotten, but it is pointed out 
that such gas component plays an important role in ex- 
plosive eruption.^ Type of volcanic eruption is classified 
into three classes by chronological behavior; One is so- 
called Vulcanian type eruption, which is widely observed 
in Japanese volcanos. This type is characterized by an 
intermittent explosive eruption and formation of a lava 
dome. These features are determined by physical prop- 
erties of magma; Specifically viscosity of magma controls 
them. 

In this paper, we study Vulcanian eruption, because its 
explosive mechanism will be the most interesting phys- 
ically, in particular, in the context of nonequilibrium 
physics; In the volcanology, an eruption picture is consid- 
ered as follows: A stage of eruption dynamics consists of a 
magma chamber and a conduit. Top of conduit is covered 
by a lava dome. In a top of magma chamber or a lower 
part of conduit, a gas component is almost completely 
dissolved into the magma molt. In the upper region of 
saturated magma, the gas is exsolved according to the 
equilibrium solubility law. As decreasing the lithostatic 
pressure, volume fraction of gases is increasing. At the 
beginning of eruption, pressure of the magma-gas mix- 
ture is considered to increase, although the mechanism 
is not clear yet. When the lava dome cannot support 
this overpressure, it disrupts the lava dome. At the next 
moment, two shock waves appear and propagate; One 
is a shock wave formed between atmosphere and com- 
pressed air and it propagates upward. Another is decom- 
pression wave in magma-gas melt and it goes to opposite 
direction. During the eruption it is observed that the 
transition from the laminar flow of bubbly melt to the 
turbulent flow of gas-magma dispersion in the conduit. 
This transition layer determines the front of fragmenta- 
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tion wave which propagates downward. At the moment, 
viscosity of magma-gas mixture is drastically changed 
abruptly about the order of 10^^ ~ lO^^Pa • s. 

There have been many theoretical investigations of 
Vulcanian eruption in the volcanism study. In 1995, 
Woods proposed the model for magma flow in conduit;^ 
In his model magma-gas mixture is treated as a one- 
dimensional nonviscotic compressible fluid with single 
component. This model can capture physical properties 
of dynamics in some sense. But treatment of dynamics 
is not well satisfied; For example, flow is treated as isen- 
tropic one, though bubble nucleation accompanies the 
eruption. There are some other phenomenological mod- 
els, but the present understanding of the eruption dy- 
namics is still unsatisfactory in the context of nonequi- 
librium physics.^' 

Recent progress of experimental techniques enables us 
to compare such theoretical model with experimental re- 
sults; These experiments are called as shock-tube ex- 
periment.^"^ In the experiment, analogue materials of 
magma-gas mixture, such as viscoelastic materials and 
powder, are used. It is observed that the behavior of ex- 
plosion depends on the viscosity of analogue materials. 
Thus a non-viscotic treatment in a theoretical study is 
not sufficient. 

In this paper, we try to establish a computational mi- 
croscopic model of Vulcanian eruption; So to say, we want 
to make "an Ising model of Vulcanian eruption" . Here 
wc describe dynamics of the mixture by microscopic par- 
ticle dynamics. A particle dynamics simulation can be 
regarded as an ideal shock-tube experiment, because we 
can calculate macroscopic quantities. In addition, using 
the particle dynamics, we can also reproduce hydrody- 
namic behavior described by a continuum description of 
Navier-Stokes equation. Even in Newtonian dynamics, 
we can produce macroscopic behavior in linear nonequi- 
librium thermodynamic regime. "'^''"^^ Moreover we can 
also discuss phenomena in far from equilibrium state, 
which are not captured by continuum descriptions based 
on local equilibrium. Thus the particle model enable us 
to explore nonequilibrium dynamics of volcano, as well as 
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Fig. 1. Geometry of the system. When we calculate physical 
quantities, we slice the system with a unit length. 

the model can verify an macroscopic theoretical model. 

Here we assume microscopic dynamics are governed by 
the following Hamiltonian: 

N 2 1 ^ 

1=1 ' i,j 

where 4>{r) is Lennard- Jones 12-6 potential: (/)(r) = 
4e{(CT/r)^^ — (cr/r)^} + 0o- For computational efficiency, 
we introduce a potential cutoff as 3.9(t and determine 
the value of </io to be (/)(3.9o') = 0. And N denotes to- 
tal particle number, rrii denotes mass of particle i, 
and (\i denote particle three-dimensional momenta and 
coordinates, respectively. Dimcnsionless parameters Ui 
and ?Tigas/?Tiniagma are selected so that it will reproduce 
similar properties as magma gas;^ We take a; to be 1 
for magma particles, and 0.1 for gas particles. It deter- 
mines energy scales of magma and gas. Ratio of melting 
temperatures of magma to gas is given by ctmagma/'^gas 
and it is 100 in the present model, although it is ap- 
proximately 1000 for actual magma and gas. Present 
choice is ten times less than actual situation, but it is 
sufficient to describe the explosive eruption as we will 
show in the following. Particle mass ratio is chosen as 
'7^gas/™magma = 0.1, which is of Order actual mass ra- 
tio. Hereafter we measure length, mass, and energy by 
the units of cr, mmagma and e, respectively, and use di- 
mcnsionless variables. Employing the Lennard- Jones 12- 
6 potential makes us to describe thermodynamic phases 
of gas, fluid, solid, and their coexisting state. 

Using the above Hamiltonian, we calculate particle 
motion. The geometry of the system is as follows (see also 
Fig. 1): Consider rectangular parallelepiped with a size 
Lx X Ly X Lz- For x and y directions, periodic boundary 
conditions are imposed. A eruption direction is to z axis, 
and we prepare elastic walls at bottom and top. These 
walls are represented by repulsion part of Lennard-Jones 
potential. 

First we have to prepare initial state as thermal equi- 
librium one. In this stage, whole system is divided into 
two parts, "chamber" {0 < z < L^) and "conduit" 
[Ld < z < Lz) hy a, diaphragm, which is located at 
z = Ld, made of same elastic walls at z = and z = L^- 
At the beginning, magma and gas particles are contained 
in the chamber. Contrarily, only gas particles are in the 
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Fig. 2. Space-time profile of number density (left) and local pres- 
sure (right): Horizontal axis represents coordinate of explosion 
direction {z axis) and vertical axis is time. At time 0, a di- 
aphragm is removed. Characteristic waves are guided by lines. 

conduit. For preparing initial state, we do an isothermal 
simulation with Nose-Hoover thermostat in each part of 
the system. ^"^"^^ Density and temperature in the cham- 
ber are chosen as gas particles are uniformly mixed into 
magma particles; There is no phase separation. 

After thermalization, we remove the separator between 
conduit and chamber and we detach the thermostat. 
Then the system obeys the Hamiltonian dynamics. If 
pressure in the chamber is higher than one in the conduit, 
an explosion is activated. 

Simulation details are as follows: The second order 
symplectic method (the leapfrog method) is used in nu- 
merical integration. Time integration slice is taken to 
be 10~'^. This value is sufficient for present simulations, 
which is checked by energy conservation. 

In the simulation, we calculate several physical quan- 
tities in boxes which are obtained by slicing along z- 
direction with a unit length cr.^''' Number density n{z) 
and mass density p(z) of the slice z are basic quantities 
of macroscopic dynamics defined by counting a number 
and mass in the local slice. Barycentric velocity v(z) is 
defined through sum of momenta in the slice. Pressure 
p{z), is defined by a trace of stress tensor. And temper- 
ature T{z) is defined by variance of particle velocities 
from local barycentric motion. 

Here we present a typical result of simulation as space- 
time profile of physical quantities. In Figs. 2, number 
density n{z) and pressure p{z) are presented. In this 
simulation, we take following parameters: System size is 
Lx = Ly = AO, Lz = 740. In an initial thermal equilibra- 
tion stage, a diaphragm is located at z = 40, so the size 
of magma chamber is 40 x 40 x 40 and one of the con- 
duit is 40 X 40 X 700. Total number of particle is 176 000, 
which consists of 57 600 magma particles and 118 400 gas 
particles. The chamber contains 57 600 magma particles 
and 6 400 gas particles. Other 112 000 gas particles are 
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Fig. 3. (Color online) Snapshots of simulation: (Up) Snapshot at 
t = 40. (Down) Snapshot at t = 170. Parameters are identical to 
ones of Fig. 2. Eruption propagates to the right direction. Only 
particles originated from the chamber are plotted; A red ball 
represents a magma particle, and blue one is a gas particle. At 
the initial condition t = 0, blue and red particles are uniformly 
mixed in the chamber. 

in the conduit. Then initial number densities are 1 for 
the chamber and 0.1 for the conduit. Thcrmalization is 
done with the chamber temperature 2 and the conduit 
temperature 0.8. 

In Figs. 2, a horizontal axis corresponds to z direction 
and explosion goes to right. A vertical axis represents 
time. At the time 0, the diaphragm is removed. In the 
profile of number density, we recognize two characteris- 
tic density waves. First one begins at {z = 40, t ~ 0) 
and propagates to (750, 120). This wave corresponds to 
a shock wave between hot gas, which is heated by adia- 
batic compressing, and thermal equilibrium gas. Its ve- 
locity is larger than a sound velocity of equilibrium con- 
duit gas. This wave is reflected at (757,120), because 
an elastic wall exists at there. Another wave propa- 
gates more slowly than the shock wave from (40,0) to 
(300,185). Front position of this density wave corre- 
sponds to magma-gas contact surface. 

There are other small waves in this figure. A wave 
propagating from (0, 10) to (170, 185) is also reflecting 
wave caused by the elastic wall located at z = 0. A wave 
propagating to opposite direction, which is from (40, 0) 
to (0, 10), is also observed in the figure. This wave is an 
expansion wave of dense magma-gas mixture. 

Other significant features are observed in this space- 
time profile. After propagating magma-gas contact wave, 
some internal structures are glowing. To investigate the 
internal structure in details, we show snapshot of sim- 
ulation are shown in Figs. 3. These snapshot are taken 
from the simulation drawing Fig. 2, so simulation pa- 
rameters are identical ones of that simulation. We only 
draw magma particles and gas particles which are in the 
magma chamber at the initial condition. Gas particles 
coming from the conduit arc omitted. Explosion propa- 
gates to the right direction in this figure, which is z axis. 

Before removing the diaphragm, magma and gas are 
uniformly mixed in the magma chamber. But, in Figs. 3, 
inhomogeneous mixing of those components is gradually 
growing during the eruption. This reminds us of spinodal 
decomposition. Size of exsolved gas bubble grows from 
Figs. 3(a) to (b); In Fig. 3(a), bubble size are widely 
distributed but, in (b), one large gas bubble and small 
bubbles in the thick magma exists. In large gas bubble. 



one magma droplet is observed. 

In this way, magma-gas mixture become inhomoge- 
neous mixture and internal structure of bubbles are grow- 
ing. Such behavior is consistent with the scenario of vol- 
canology. But in the present simulation, transition to 
magma dispersion fiow is not observed. The reason may 
be that smaller cross section of conduit and finitcness 
particles. 

Next we compare the present simulation results 
with the continuum description given by Woods. ^ In 
his model, magma-gas mixture is described by one- 
dimensional nonviscotic compressible one-component 
fluid. The dynamics are described by a continuity equa- 
tion, an equation of motion, and the foUowings: 

1-71 nRT 1 f (jiV"^ 

\ = -, Pg\-] const. , (2) 

Pi Pg P \PJ 

where p, pi,pg,T, R,n, (j) and denote mass density, 
mass density of magma component, pressure of gas com- 
ponent, temperature, a gas constant, a mass fraction 
of magma and gas components, a volume fraction of 
magma and gas components, and ratio of specific heats, 
respectively. In these quantities, p,pg,T and (/)~^ = 
1 -I- ■^-^ p^^rp are variables. Other pi,R,n, and 7,„ are 
fixed to some constant values. The first equation is an 
equation of states, and the second one expresses an isen- 
tropic condition derived from the first law of thermo- 
dynamics. As we know the present equation of states is 
almost identical to one of ideal gas. 

These equations are essentially same as ones of com- 
pressible ideal gas fluid. To study the equations is just 
an textbook example. ^^'^^ We get a standard rewrite as 

where w and a{p) are a velocity field and a sound veloc- 
ity, respectively. The sound velocity of magma-gas mix- 
ture is a function of p, and it is expressed as a^(p) — 
al{p/poV"'''^{(l)o/(f'V"'^^^ (ao,Po,^o are sound velocity, 
density, and volume fraction at some reference state.) 
This equation gives characteristic curves and conserved 
quantities on them. Then we can solve the equation in 
characteristic regions. For obtaining global shock tube 
solution, we have to glue the solution with appropriate 
boundary conditions. 

In Fig. 4, temperature T{z), barycentric velocity 
v(z)^, pressure p{z), mass density p{z) of the present 
simulation are shown. Simulation parameters are taken 
to be as follows: System size is — Ly = 32, — 408 
and size of magma chamber is 32 x 32 x 200. Initial num- 
ber density of magma chamber is taken to be 1 and con- 
duit density is 0.02, thus the number of particles in the 
chamber is 204 800, which contains 10% gas particles. 
The number of gas particles in the conduit is 4 096. In 
this simulation, we imposed an artificial boundary condi- 
tion at the top of conduit; For decreasing reflection effects 
from the top elastic wall, we attach a particle sink at the 
top, in which particles with the energy larger than some 
threshold value are removed from the system. 

We can observe characteristic regions in Fig. 4. Let us 
compare these results with continuum descriptions. The 
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Fig. 4. Spatial profiles of temperature, velocity (z) , pressure, and 
mass density at t = 15: System size is taJjen to be Lx = 32, Ly = 
32, Lz = 408 and size of magma chamber is 32 x 32 x 200. Ini- 
tial mass density and temperature are taken to be 1 and 2, re- 
spectively. We can recognize characteristic regions. From right, 
"initial equilibrium state" , "hot gas region" , "cold gas region" , 
"expanding wave region", and "initial equilibrium state" again 
are observed. These regions are indicated by gray rectangular. 



solution of Eq. (3) teaches us that there are three regions 

in shock tube analysis, that is, a hot gas region, a cold 
gas region, and an expanding wave region. Corresponding 
regions of simulation are indicated in the figure; In the 
"hot gas" region, gases are heating up by the shock wave. 
In contrast, in the "cold gas" region, gases are cooling by 
an adiabatic expansion. Another region is an "expand- 
ing wave" region in which the expanding wave exists and 
physical quantities are smoothly changed. Physical prop- 
erties of such regions obtained by the simulation are al- 
most equivalent to ones of shock tube analysis. But there 
is a little mismatch with the solution; Analysis of com- 
pressible fluid gives constant profiles of physical quanti- 
ties in both hot and cold gas regions. But, in this simu- 
lation, some structures arc observed in each regions. For 
example, in a velocity profile of the hot gas region, veloc- 
ity near cold gas is rather faster than other areas. This 
high velocity area is caused by pushing effects of magma- 
gas contact surface, which is corresponding to the front 
of cold gas contact. These high velocity particles are not 
thermalized yet; In the molecular dynamics simulation, 
microscopic relaxation is apparently observed. 

To summarize, we have constructed a microscopic 
model of Vulcanian eruption by a two-components 
Lennard- Jones particle system. We observed that the 
particle dynamics is efficient in this kind of dynamics. 
Using the present model we can reproduce characteris- 



tic features of explosive eruption such as a shock wave, 
a expansion wave. At the early stage of the eruption, 
we also compare the simulation result with the analytic 
model given by Woods. Qualitative behavior is almost 
consistent with the analytic result, even though the flow 
is treated as nonviscotic one in the analytic model. In ad- 
dition, we have also observed that the internal structure 
is growing during the eruption. Internal bubble structure 
cannot be captured by the Woods model. This behavior 
is also consistent with a eruption picture of volcanology 
study. Thus we conclude that the present model is a can- 
didate of "an Ising model of Vulcanian eruption" . 

To establish the present model, a quantitative study 
is inevitable. For this purpose, we have to enlarge the 
size of system; A transition from bubbly magma flow 
to magma dispersion flow will be reproduced and stud- 
ied by simulation of the system with ten times larger to 
all directions. And the more details of volcanic eruption 
not only Vulcanian but also Strombolian, and Plinian 
will be elucidated. Present typical computational time is 
approximately 80 hours for 208 896 particles with single 
AMD opteron 248 (2.2GHz). Hence much larger simula- 
tion is feasible with large super computers. 
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